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Abstract. One of the key questions to understanding 
the efficiency of diffusive shock acceleration of the cosmic 
rays (CRs) is the injection process from thermal parti- 
cles. A self-consistent injection model based on the inter- 
actions of the suprathermal particles with self-generated 
magneto-hydrodynamic waves has been developed re- 
cently by Malkov (199J). By adopting this analytic so- 
lution, a numerical treatment of the plasma-physical in- 
jection model at a strong quasi-parallel shock has been 
devised and incorporated into the combined gas dynam- 
ics and the CR diffusion-convection code. In order to in- 
vestigate self-consistently the injection and acceleration 
efficiencies, we have applied this code to the CR modi- 
fied shocks of both high and low Mach numbers (Af = 30 
and M = 2.24) with a Bohm type diffusion model. Both 
simulations have been carried out until the maximum mo- 
mentum (pmax/mpc) ~ 1 is achieved to illustrate early 
evolution of a Bohm type diffusion. We find the injection 
process is self-regulated in such a way that the injection 
rate reaches and stays at a nearly stable value after quick 
initial adjustment. For both shocks about lO"'^ of the in- 
coming thermal particles are injected into the CRs. For the 
weak shock, the shock has reached a steady state within 
our integration time and ^ 10% of the total available 
shock energy is transfered into the CR energy density. The 
strong shock has achieved a higher acceleration efficiency 
of ~ 20% by the end of our simulation, but has not yet 
reached a steady-state. With such efficiencies shocks do 
not become CR-dominated or smoothed completely dur- 
ing the early stages when the particles are only mildly rel- 
ativistic. Later, as the CR pressure becomes dominated by 
highly relativistic particles that situation should change, 
but is difficult to compute, since the maximum CR mo- 
mentum increases approximately linearly with time for 
this model. In the near future we intend to extend such 
shock simulations as these to include much higher CR mo- 



menta using an adaptive mesh refinement technique cur- 
rently under development. 
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1. Introduction 

The non-thermal energy distributions of cosmic ray ions 
or source distributions of electrons emitting synchrotron 
radiation in various astrophysical objects are commonly 
described as produced by the first order Fermi acceleration 
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process at shocks (for reviews see Drury 1983; Blandford 
& Eichler |198^; Kirk et aL [l994| ). 

When particles diffusetl off the moving scattering cen- 
ters in a region divided by a velocity discontinuity (shock), 
these particles can be accelerated if their mean free paths 
exceed the shock thickness. The relative momentum gain 
for a cycle of two crossings of the shock is then propor- 
tional to the velocity difference across the shock, i.e. of 
first order with respect to the shock velocity (Bell 1978). 
In astrophysical coUisionless plasmas an electro-magnetic 
field must be present to change the energy of particles. 
Waves or irregularities in this field provide particle scatter- 
ing, which leads to diffusion. Consider a shock with veloc- 
ity Us > propagating into a plasma at rest with density 
p and with a homogeneous magnetic field -Bo in the direc- 
tion of the shock normal. The plasma is compressed to the 
density pd by the shock, and flows downstream with the 
velocity Ud = Ms(1 — l/''); where r = Pd/p is the compres- 
sion ratio. Particles with the mean downstream velocity 
(v) — Ud cannot cross the shock from downstream to up- 
stream, because (v) < itg. In addition the shock may not 
be a discontinuity for a particle at this energy, because 
the gyro radius of the thermal particles is of the order 



^ In general, also anomalous transport like siifc-diffusion or 
super-diffusion can be realized. 
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of the shock thickness, leading to adiabatic energy change 
while crossing the shock. Because the plasma is also heated 
downstream of the shock, some supra-thermal particles in 
the high energy tail of the Maxwellian velocity distribution 
may gain energies and have velocities that allow them to 
re-cross the shock. In a homogeneous magnetic field, due 
to lack of scattering centers these particles would escape 
upstream without returning to the shock. Then, the accel- 
eration mechanism would not apply. However, the popu- 
lation of particles that can move upstream provide a seed 
particle beam, which generates Alfven waves responsible 
for scattering and, therefore, diffusion, which is an essen- 
tial element of the first order Fermi acceleration process. 

The problem of particle acceleration from thermal en- 
ergies up to relativistic particle energies is highly non- 



linear, as first pointed out by Eichler (1979). First, the 
energy transferred from the bulk of the plasma to the 
sub-population of accelerated particles can change the 
thermodynamic properties of the plasma like the temper- 
ature and density. In addition, the accelerated particles 
provide their own pressure in the system, which, since 
it differs from the thermal pressure, modifies the veloc- 
ity structure of the shock transition. Second, the waves 
generated by particles escaping upstream determine the 
transport properties of the plasma, and, therefore, reg- 
ulate this wave generating escape itself. The manner in 
which the wave-particle interactions control the fraction 
of plasma particles that can escape upstream to partic- 
ipate in the Fermi process is commonly called injection. 
This is a basic aspect of the plasma of coUisionless shocks 
and is itself highly non-linear. This injection problem is 
fundamentally related to the question of the efficiency of 
particle acceleration at shocks by the Fermi process. 

Different numerical methods have been used to treat 
the injection problem of CR modified shocks. In Monte- 
Carlo simulations of non-linear particle acceleration the 
details of a posited scattering law provide an injection pa- 
rameter, but one not determined self-consiste ntly f rom the 
partic le wa ve interaction (e.g. Ellison et al. 1996 ; Baring 
et al. 1999 ). In contrast to pure kinematical effects from 
shock velocity, particle speed and inclination angle of mag- 
netic field and shock, the waves responsible for particle 
scattering depend on plasma properties like temperature 
and the beam strength of the wave generating particles 
itself. These coupled and time dependent effects are not 
easy to incorporate into a Monte-Carlo approach. Time 
dependent Monte-Carlo simulations have been presented 
by Knerr et al. ( 1996| ), but still with a prescribed scatter- 
ing law, as a parameterization of injection. 

In the two-fluid approach the cosmic rays are treated 
as a diffusive gas without following their momentum dis- 
tribution. The energy transfer into CRs in these models 
is based on a fraction of the upstream gas particle s, tha t 
are instantaneously accelerated a t the shock (Dorfi |l990 ), 
around the shock ( Jones & Kang 1990 ) or at velocity gra- 
dients (Zank et al. 1993 ). In practical terms, because the 



shock is the most prominent velocity gradient in the sys- 
tem, these techniques are very closely related, as pointed 
out by Kang & Jones ( 19951 ). Essentially the same param- 
eterization is also used in the numerical solution of the 
hydro-dynamical equations coupled to the momentum de- 
pende nt co smic-ray trans port e quation (e.g. Falle & Gid- 



dings |l987t Kang & Jones |1991| ; Berezhko et al. |l994| ). 

Kang & Jones (1995) used a numerical injection 
model with two essentially free parameters which describe 
boundaries in momentum at which particles can be accel- 
erated {pi — ci -pth where Pth is the peak momentum of the 
Maxwell distribution) and from which these contribute to 
the cosmic-ray pressure (p2 ~ C2 -Pth). Still, these momen- 
tum boundaries are free parameters, which can be trans- 
lated into a particle fraction of the upstream gas. How- 
ever, models that incorporate more details of the plasma 
physics of the background plasma, and, therefore, the CRs 
at injection energies are really necessary to constrain the 
phase-space function and therefore determine these pa- 
rameters. In this way we can incorporate self-consistent 
plasma physical models into numerical simulations of CR 
acceleration. 

Such a plasma physical model based on non-linear 
interactions of particles with self-generated waves in a 
shocked plasma has been investigated numerically by solv- 
ing the kinetic equations of ions in a magnetic field, and 



treating the electrons as a background fluid (Quest 1988) 



These simulations show that ions can be scattered back 
and fo rth ac ross the shock by self-generated waves, and 
Quest ( 198^ also points out that these scattered ions can 
provide a seed population of cosmic rays. 

Recently, the kinetic equations of ions were solved an- 
alytically for non-linear wave-fields near strong paral- 



lel shocks by Malkov & Volk ( |1995D and Malkov (|998|). 
These authors were able to constrain the fraction of phase- 
space of the background plasma that can be injected into 
the acceleration process as a result of the self-regulating 
interaction between wave generation and particle stream- 
ing. Here we incorporate this self-consistent analytical re- 
sult in numerical solutions of the hydro-dynamical equa- 
tions together with the cosmic-ray transport equation. 
Our simulations, therefore, provide the first time depen- 
dent solution of the problem of modified shocks that in- 
cludes a self-consistent plasma physical injection model. 
This technique enables us to determine the level of shock 
modification and acceleration efficiency in an evolving 
shock without a free parameter for the injection process 
(Gieseler et al. |1999D . 

We describe the plasma physical injection model in 
some more detail in Sect. ^ before we outline the coupled 
set of dynamical equations of the plasma and cosmic rays 
in Sect. |[ together with details of the numerical method 
we used to solve these equations. The results are presented 
in Sect. || and Sect. |[ consisting of the time dependent 
evolution of the plasma properties, showing especially the 
modification of the velocity profile and the momentum dis- 
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tributions of thermal and rclativistic cosmic-ray particles. 
In addition we present our results in the form of a parti- 
cle injection efficiency and an energy transfer efficiency at 
modified strong shocks. 

2. Injection Model 

2.1. Wave-particle interaction at parallel shocks 

In supernova remnants (SNRs) the relative orientation of 
magnetic field and shock front can be very diverse, even in 
one single object. For example, if a spherically symmetric 
shock front expands in a region of homogeneous magnetic 
field, the directions of the shock normal and the magnetic 
field change over the shock surface from parallel to perpen- 
dicular. For nearly perpendicular shocks the acceleration 
process can be very fast and effective due to ref lections up- 
stream of the shock (Naito & Takahara |1995| ). However, 
the velocity of the intersection point of shock and mag- 
netic field in highly oblique shocks can be close to light 
velocity. This can suppress the injection efficiency of ther- 
mal particles, which are effectively tied to magnetic field 
lines. That is because the velocity distribution function of 
thermal protons, which we shall assume to be Maxwellian, 
drops sharply towards high energies. This purely kinemat- 
ical effect has been investigated by Baring et al. (1993). 
Therefore, regions of SNRs where quasi-parallel shocks 
exist, are likely to be where the most effective injection 
occurs. On the other hand, effective acceleration, i.e. short 
acceleration time scales and hard spectra, may be realized 
in other parts of a SNR, where an oblique geometry of 
magnetic field and shock normal is found. 

For quasi-parallel shocks, where the shock propagates 
along the mean magnetic field (x-direction) , the transport 
properties along the mean field direction are most impor- 
tant. We will assume this case, with the field Bq paral- 
lel to the shock normal. The spatial diffusion of particles 
is produced by magneto-hydrodynamic waves, which are 
in turn generated by particles stre aming along the mag- 
netic field, Bq. We refer to Malkov ( 1998 ) for an extended 
analytical description of the particle-wave interaction for 
low-momentum particles, and we describe here only the 
results which are relevant for the implication of this model 
in our simulations of the time dependent acceleration at 
modified shocks. When particles are streaming along the 
magnetic field in the upstream direction, waves are gen- 
erated due to the ion cyclotron instability. The resulting 
upstream magnetic field, which corresponds to a circularly 
polarized wave, can be written as 



B = BqBx + [By cos kox — Bz sin kox) 



(1) 



The amplitude will be amplified downstream of the 
shocks by a factor r = p^/ p. The downstream field can be 
described by a parameter e, for which, following Malkov, 
we assume e := Bq/B^^ ^ 1, in the case of strong shocks. 
Note that the perpendicular component of the magnetic 




downstream 



upstream 



Fig. 1. Cartoon of the injection model in the shock-frame 
phase-space. Plasma is moving towards —x into the shock with 
velocity —its and gets compressed, heated and decelerated to 
the downstream velocity —U2. Particles with positive veloc- 
ity can stream back to upstream along the magnetic field Bq. 
These particles provide the beam, which generates the mag- 
netic field waves. The magnetic field wave is shown schemat- 
ically in configuration space. The wave amplitude, frequency, 
and damping length is shown only qualitatively. 



field leads effectively to an alternating field downstream 
of the shock for particles moving along the shock normal 
(see Fig. I). 

2.2. Thermal leakage model 

The particles with a large enough gyro radius 



r^ 



— pc sinQ;/(ei?^) > I/Zcq 



(2) 



can have an effective velocity with respect to the wave 
frame, i.e. the downstream plasma would be transparent. 
Some of these particles that are in the appropriate part 
of the phase space (depending on the shock speed) would 
be able to cross the shock from downstream to upstream. 
For the protons of the plasma, the resonance condition 
for the cyclotron generation of the Alfven waves gives 
kQ{v) « = u!±Bq/B^, where the cyclotron frequency of 
protons is given by uj^ = eBj^/{mpc), and {v) is the mean 
downstream thermal velocity of the protons. We now have 
for the thermal protons korg± ~ e ^ 1. This means that 
most of the downstream thermal protons would be confined 
by the wave, and only particles with higher velocity in the 
tail of the Maxwellian distribution are able to leak through 
the shock. Ions with mass-to-charge ratio higher than pro- 
tons have a proportionally larger gyro radius, so that the 
injection efficiency of protons would yield a lower limit 
for the less magnetized ions. On the other hand, for ther- 
mal electrons a plasma with such proton generated waves 
would have a reduced transparency due to the smaller 
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Fig. 2. Transparency function Eq. (3) vs. the normalized par- 
ticle velocity v = vko/u!± for e = 0.35 (solid line). Shown as 
a dashed line is the transparency function given in Fig. 2 of 
Malkov & Volk (1998). 




Fig. 3. Transparency function Eq. (3) vs. the normalized par- 
ticle velocity v = v/u2 for different values of e. We used the 
relation uj±/k(, ~ U2/e (see text). 



gyro radius of the electrons. However, reflection of elec- 
trons (and protons) off the shock could become efficient 
with increasing wave amplitude and possibly aid in their 



injection (e.g. Levinson 1996; McClements et al. 1997). In 
the following we will focus on the protons, which carry 
most of the energy and momentum of the plasma. 

2. 3. Transparency function 

To find the part of the thermal distribution for which 
the magnetized plasma is transparent, a nd, wh ich, there- 
fore, forms the "injection pool", Malkov ( 1998 ) solves an- 
alytically the equations of motion for protons in self- 
generated waves. He finds a transparency Junction t^^^, 
which expresses the fraction of particles that are able to 
leak through the magnetic waves, divided by the part 
of the phase space for which particles would be able to 
cross from downstream to upstream when no waves are 
present. For the adiabatic wave particle interaction the 
transparency function is given by Malkov (1998) Eq. (33), 
with r^gc = 21^^30/(1 ~ U2/v), where v is the particle ve- 
locity V and U2 = Ug/r is the velocity of the shock in the 
downstream plasma frame. Here v^^^ is the fraction of the 
particles streaming back from downstream to upstream. 
This quantity is divided by the fraction of particles that 
would be able to escape upstream in the absence of waves. 
In order not to further increase the complexity of our nu- 
merical simulation, we use here the following approxima- 



tion of the representation given in Malkov (1998) 



tUv, U2) = H[d-il + e)] (1 - ^) ' (1 - ^ 

:{-[v-il + e)]-'] , (3) 



• exp ■ 



where the particle velocity is normalized to v = vko/cu^ 
and H is the Heaviside step f unctio n. We argued above 
that Lu^/ko ~ u2/e (see Malkov |l998| Eq. 42). The trans- 
parency function now solely depends on the shock velocity 



in the downstream flow frame, U2, the particle velocity, w, 
and the relative amplitude of the wave, e. 

The calculation of the transparency function and the 
wave-amplitude e uses the ergodicity of the downstream 
phase-space for the randomized motion of particles in 
the high-amplitude wave field. The upstream wave field 
is generated by a beam of leaking particles whose en- 
ergy density is calculated from the corresponding area 
of the downstream phase-space. From the energy den- 
sity of this beam the upstream magnetic-field wave ampli- 
tude is determined self-consistently (Malkov 1998 ).□ Be- 
cause of this feedback Malkov was able to constrain the 
quantity e as 0.3 ^ e ^ 0.4, leaving essentially no free 
parameter. Comparison with hybrid plasma simulations 
suggests 0.25 < e < 0.35 (Malkov & Volk [l998[ ), con- 
sistent with their analytical results. This constraint on 
the wave-amplitude B^/Bo = 1/e defines the level to 
which the particle-wave interaction adjusts. With the es- 
timation of this amplitude there is no free parameter de- 
scribing the level of the beam strength for the injection, 
and, therefore, the injection efficiency. The advantage of 
our approach presented here is that quantities like the 
plasma velocity and particle momentum distribution are 
calculated self-consistently by solving simultaneously the 
hydro-dynamical equations together with the cosmic-ray 



transport equation (see Sect. 3.1). 

The function (3) is plotted in Fig. || for e — 0.35 vs. the 
particle velocity normalized to v = vko/u!±. In Fig. |^ we 
have also reproduced this function as given in Fig. 2 of 
Malkov & Volk ( |1998D to allow a direct comparison. The 



^ Especially downstream, where Bo/B± ^ 1, the particle- 
wave interaction can by no means described by quasi-linear 
theory, which would be valid in the opposite case, i.e. for in- 
coherent waves with small amplitudes. In fact, the calcula- 
tion of the transparency functio n is n ot based on results of 
the quasi-linear theory (Malkov 1998| ; Malkov, private com- 
munication) . This point is noteworthy, becau se of existing mis- 
interpretations of the work of Malkov (199J) in the literature 
(Baring [1999I). 
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strong velocity dependence and also the asymptotic be- 
havior is modeled reasonably well by the representation 
Eq. (3). In the normalization of Fig. ^ the dependence on 
e is very weak, and, therefore, not shown. To illustrate the 
dependence of the transparency function on small varia- 
tions of the field amplitude, it is better to choose a differ- 
ent normalization. Therefore, the transparency function 
Eq. (3) is shown in Fig. ^ vs. the velocity normalized to 
V — v/u2 for the maximal allowed range in e, as described 
above. 

3. Model 

3.1. Dynamical equations 

The standard hydro-dynamical equations of mass, mo- 
mentum and energy conservation for a gas with velocity 
u{x,t), and density p{x,t), corrected for CR pressure ef- 
fects are given by 



dp 

dt 
Au 



P 



dt 
dt 



du 
^ dx ' 



d_ 

dx 



Pc), 
u(Pg+Pe)] 



P.|-5(.,t), 



(4) 
(5) 
(6) 



where Pg and Pc are the gas and the CR pressure, respec- 
tively, and eg = Pg/ p(7g — l)-|-u^/2 is the total energy den- 
sity of the gas per unit mass. Here d/dt = d/dt+ud/dx is 
the total Lagrangian time derivative. We assume 7g = 5/3 
for the thermal gas adiabatic index throughout this work. 
The injection energy loss term S{x, t) accounts for the en- 
ergy transferred to high energy particles and will be dis- 
cussed later. Equations (4)-(6) are solved using a Total 
Variation Dimin ishing (TVD) code based on the scheme 
of Harten (|1983|) . 

We assume that the shock Mach number M = Wg/cs 
(with Cs = (pjPlp)^!'^ the upstream sound speed) ex- 
ceeds the Alfven Mach number A/a — Us/ca <C M (with 
Ca = i3/(47rp)^/^ the upstream Alfven speed). Then the 
diffusion-convection equation, which describes the time 
evolution of the phase-sp ace d ensity f{p,x,t) of the high 
energy CRs (e.g. Skilling 1975), takes the form: 



df Idu df d f , , a , 
dt 3dx dp dx \ ' dx 



(7) 



The diffusion coefficient k{p, x) is assumed to be a 
scalar. Transforming to the variables y := ln(p ) and 
g{y,x,t) -.^ p ' ^ f{p,x,t) {cf. Falle & Giddings |l987| ; Kang 
& Jones [l99l] ) , Eq. (7) can be written as 



dg 
dt 



1 du 
3 dx 



dg 
dy 



45 



d_ 

dx 



(8) 



This equation is solved using an implicit Crank-Nicholson 
scheme, which is sec ond order in space and time (see 
e.g. Falle & Giddings |l98l. 



The high energy particles provide an additional pres- 
sure to the system that has to be included in the set of 
hydro-dynamical equations (with p normalized to the pro- 
ton momentum p/vripC —>■ p): 



(9) 



This definition of the CR pressure Pc includes the dif- 
ference of the phase-space density from the Maxwellian 
distribution and defines the sub-population which we 
identify as CRs. The CR energy density is defined as 



= Anm^c' J (/ - /m)p'(Vp' + 1 - l)dp. 



(10) 



3.2. Injection scheme 

We do not include an additional injection term in 
Eq. (8), because in our model injection is described self- 
consistently from the thermal distribution. Therefore, the 
lower boundary of the momentum distribution of the CR 
population must match the upper boundary of the mo- 
mentum distribution of the gas. The distinction between 
these populations is, of course, only technical, and defined 
by the validity of the relevant dynamical equations. We 
use a Maxwell distribution according to the actual den- 
sity and temperature of the plasma. Instead of a fixed 
momentum boundary we use here the transparency func- 
tion T„s^ to define where the lower boundary of the CR 
momentum distribution matches the momentum distribu- 
tion of the bulk plasma. The injection into the high energy 
part of the phase-space distribution (i.e. the CRs) is then 
directly provided by the bulk of the plasma. 

The initial Maxwellian phase-space density fmip, x, t) 
is given by: 



i{p,x,t) = p'^f^.i{p,x,t) 
n(x, t)p'^ 



(27rmp k^Tf/^ 



exp 



2?7ip k^T 



(11) 



where n{x,t) — p/rup is the particle number density, and 
the temperature is defined by the local gas pressure Pg 
and density p according to T = prUpPg/ pk^. Here p is 
the mean molecular weight which is assumed to be one, 
and fee is the Boltzmann constant. The details of how 
the momentum distribution is calculated in a time step 
from t — At to t are as follows. First we define the CR 
part of the momentum distribution by gcnit — At) = 
g{t — At) — gm{t — At). Now the CR diffusion-convection 
equation (Eq. 8) is solved for the entire momentum space, 
including the thermal Maxwell distribution, to find the 
updated distribution function g(j),x,t). For momenta be- 
low the critical momentum of T^sc{Pcnt) = any particle 
acceleration must be suppressed, and therefore the result 
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of Eq. (8) is rejected by restoring the Maxwellian distri- 
bution given in Eq. (11). For momenta above the criti- 
cal momentum and upstream of the shock, we use the 
transparency function as a filter for g{p, x, t) as described 
below, since t^^^ {p, t) corresponds to the fraction of the 
phase-space density at a given momentum that can cross 
the shock from downstream to upstream. The final distri- 
bution at time t is then given by g{p,x, t) in the following 
way: 

gcR{p,x,t) = g{p,x,t) - g^{p,x,t) , (12) 
gcR{p,x,t) = gcR{p,x,t - At) + Tcsc{p,t) 

■ [9cr{p, X, t) - gcn ip, X, t - At)] , (13) 
g{p,x,t) = gcR{p,x,t) + gM{p,x,t) . (14) 

So effectively at the lower momentum limit where t^^^ (p) = 
the Eq. (8) has no effect at all, while at the higher 
momentum limit where r^^cip) = 1 the result of Eq. (8) 
is used without further modification. Only in the inter- 
mediate momentum regime where < T„g^(p) < 1, the 
transparency function represents the injection process (i.e. 
thermal leakage). Thus the transparency function defines 
self-consistently the momentum boundary above which 
the particle acceleration mechanism can work, and also 
defines the transition region between thermal plasma and 
accelerated particles. 

The particle injection rate into the CR population can 
be estimated from the adiabatic change of the momentum 
due to the velocity gradient of the flow: 



f dp' 

Q{p,x,t) ^ 4TTp'^f{p,x,t) ( — 
47r o „ , , 9u 

= --p^np,x,t)-. 



(15) 



Then the energy loss rate of the gas can be written as 

oo 

Six, t) = impc^ j ^I^^^p^Qip, X, t) dp , (16) 



oo 


Note here the condition dT^^^{p,t)/dp ^ in fact de- 
fines the "injection pool" where the thermal leakage takes 
place. Due to the steep dependence of both the Maxwell 
distribution and the transparency function on the par- 
ticle momentum, the momentum range of the injection 
pool is well restricted. Either below or above this mo- 
mentum range t,^c(p) = constant, so dT^^^{p,t)/dp = 0. 
If the transparency function is given by a step function 
Tcacip) = H{p — Pi„j), it becomes the injection scheme 
adopted by Kang & Jones ( 1995 , 1997 ) in which the in- 
jection takes place at a single injection momentum rather 
than an extended momentum range. 



The transparency function t^^^ given by Eq. (3) de- 
pends on the downstream plasma velocity, which is aver- 
aged over the diffusion length of the particles with mo- 
mentum at the injection threshold. This dependence is 
also important for the injection efficiency, and leads to 
a regulation mechanism similar to the above beam wave 
interaction. If the initial injection is so strong that a sig- 
nificant amount of energy is transferred from the gas to 
high energy particles, the downstream plasma cools, and, 
in addition, the downstream bulk velocity decreases in the 
shock frame due to the shock modification of the cosmic- 
ray population. Because the injection pool is in the high 
energy tail of the Maxwellian distribution of the gas, the 
cooling decreases significantly the injection rate. However, 
the deceleration, in turn, allows for a modest increase of 
the phase-space of particles that can be injected. This is 
expressed by the U2 dependence of Eq. (3) . This velocity 
dependence balances partly the reduction of injection due 
to the cooling of the plasma. Remarkably, these two effects 
lead to a very weak dependence of the injection efficiency 
on e in the vicinity of e w 0.35. 

3.3. Diffusion model 

Since the injection process is included self-consistently, 
the diffusion coefficient is the only remaining free pa- 
rameter in our model. We assume the particle diffusion 
is based on the scattering off the self-generated waves 
which have a field component perpendicular to the plasma 
flow. The compression of the plasma leads to an am- 
plification of these waves, which is described by scal- 
ing the diffusion coefficient as k oc 1/p. In our one 
dimensional model we have to describe diffusion along 
the mean magnetic field. The lower limit for the diffu- 
sion coefficient is the Bohm diffusion coefficient Kb{p) = 
(3 • 10^2 cmVs/B^c) pV(l+p^)^/^ where B^q is the 
magnetic field strength in units of micro-gauss. For the 
present calculations we assume the diffusion coefficient is 
simply related to Bohm diffusion as 



k{p) = C Kb Po/p{x) 



(17) 



We have introduced the factor ( to account for the higher 
diffusion in the direction of the mean magnetic field, be- 
cause this direction is parallel to the shock normal, and, 
therefore, relevant for the acceleration process at quasi- 
parallel shocks. Although the time scale for the cosmic- 
ray acceleration does depend on the diffusion coefficient, 
the basic self regulation process for the injection problem 
which we investigate here is not dependent on the choice 
of (. Therefore, since our study intends to focus on the 
general time dependent behavior of this injection model, 
we do not include a completely self-consistent scattering 
model, where the diffusion coefficient is coupled to the 
spectrum of the Alfven waves. 



U.D.J. Gieseler et al.: Time dependent cosmic-ray shock acceleration with self-consistent injection 



7 



3.4- Initial and boundary conditions 

We assume there is no pre-existing CR population, so 
the initial particle distribution is purely Maxwellian with 
the local plasma temperature and density. We use open 
boundary conditions for the description of the thermal 
plasma in our simulations. In momentum space, the lower 
boundary is provided by the Maxwell distribution as dis- 
cussed above. At the highest momentum and also at the 
upstream boundary in configuration-space we use a 'free 
escape' boundary. Downstream we use a 'no diffusive flux' 
boundary, where the cosmic-ray density is always kept 
constant across the boundary. However, the grid size was 
chosen so large in the simulations we present here, that 
the cosmic-ray pressure is at all times essentially zero at 
both boundaries. 



4. Results for strong shocks 

First we consider a strong shock with an initial Mach 
number of M = 30. Unlike an ordinary hydrodynamic 
simulation, the simulation of the CR shock acceleration 
requires specification of three physical parameters, uq/c, 
po/rup, and the shock Mach number in addition to the 
diffusion coefficient. We adopted the following nominal 
physical scales for physical parameters: uq — 5000 km s^^, 
po/rrip = 0.03 cm-3, to = 4.0 • 10'' s, xq = 2.0 • lO^^ cm. 



go 



1.25 ■ 10 



"^ergcm 



We use C = 100 for 



the simulations presented here, and a magnetic field of 
B = 3p,G. The initial conditions are specified as follows: 
Pup = Pq, Uup = -tio, and Pg^up = 6.667 • lO^^Pgo in 
the upstream region, while pd — 3.987po, 1*2 = — 0.25uo, 
Pg,d = 0.75Pgo downstream. These values refiect the shock 
jump conditions in the rest-frame of the shock. 

We define the diffusion length and time at a given mo- 
mentum as Zd(p) = i^{p)Iuq and td(p) = k(j>)/uq. For a 
proper convergence, the spatial grid size should be smaller 
than the diffusion length of the injection pool particles 
{Id/xo ~ 0.04 in upstream for ~ 0.02). On the other 
hand, the spatial region of the calculation in upstream and 
downstream should be larger than the diffusion length- 
scale of the particles with the highest energies reached 
at the end of our simulation period {Id/xo 230 for 
Pmax ~ 2.5). So we used 51200 uniform grid zones for 
x/xq = [—250,250], with the shock initially at a; = and 
the grid size Ax/xq = 0.01 = 0.25 ■ Zd(Pi„j)/a:^o- We use 128 
uniform grid zones in log(p) for log(p) — [—3.0, 0.477]. We 
integrate the solutions until t — 240io which corresponds 
to td for Pmax ~ 2.5, so the CR particles became only 
mildly relativistic by the end of our simulations. 



Dynamical evolution 

Figure^ shows the normalized gas density p{x), gas pres- 
sure Pg(a;), plasma velocity u{x) and the cosmic-ray pres- 
sure Pc{x) over the spatial length x, for different times. 




Fig. 4. Gas density p/po, pressure Pg/Pgo, velocity u/uq, and 
cosmic-ray pressure Pc/Pgo, at times t = (dotted), t = 120 to 
(dashed) and t = 240 to (solid line). The shock Mach number 
is M = 30.0, e — 0.35 and = 100. The initial upstream gas 
pressure is P = 6.667 ■ 10~^Pgo. 



This shows clearly the basic features of the shock mod- 
ification by a diffusive component; that is, the adiabatic 
precursor compression and the sub-shock. The CR pres- 
sure Pc is responsible for the deceleration and compression 
of the plasma flow in the precursor region upstream of the 
sub-shock, which still remains strong. As a result, the gas 
is compressed to higher density downstream of the sub- 
shock. 

The cosmic-ray pressure immediately downstream of 
the sub-shock has not reached a steady state yet. The rea- 
son is that for a the non-thermal particles with a momen- 
tum distribution / oc p"* with s < 4, the energy density 
is an increasing function of Pmax- This applies even if the 
injection is shut down completely, like for an (5-function 
type injection in time, as shown by Drury (1983). We ex- 
pect that Pc will continue to increase after our integration 
time t = 240to, which leads to a significant modification of 
the shock structure and to the steepening of the power-law 
distribution of suprathermal particles. The simulations of 
such non-linear evolution, however, require much greater 
spatial region and grid zones and also longer integration 
time than what we could afford in our simulations. 

In real astrophysical shocks, the energy density is lim- 
ited by radiation losses in the case of electrons or more 
generally by particle escape due to the finite extent of the 
acceleration region. For the maximum energy of particles 
(Pmax ~ 2.5) achieved by t = 240 to in this simulation, 
neither effect is important, and, therefore, not included. 
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Fig. 5. Phase-space density g = p*/ vs. proton momentum 
immediately downstream of the sub-shock. Also shown is the 



transparency function Tobc 
t = (dotted), t = Wto 
and t = 240 to (solid hne) . 
and corresponding text. 



4-2. Energy distribution 



. Both functions are presented for 
(dot-dashed), t = 120 to (dashed). 
For the parameters used see Fig. 4 



The phase-space distribution g{p,x,t) = f{p,x,t) 
immediately (three zones) behind the sub-shock is shown 
in Fig. H for three different times. Initially this distribu- 
tion is given by a Maxwell distribution, as shown by the 
dotted line. At the thermal part of the distribution the 
cooling of the postshock gas due to the energy flux into the 
CR particles is responsible for the shift of the Maxwellian 
distribution towards lower energies. We have also plot- 
ted the transparency function r^^^ at the same simulation 
times. According to Eq. (16) the injection rate into the 
non-thermal distribution depends on overlap of dr^^^/dp 
and gip) that determines the injection pool. One can see 
that initially the injection rate is high and so the post- 
shock gas cools quickly, resulting in narrowing down of 
the injection pool. This causes the injection rate to de- 
crease. But then the transparency function also shifts to- 
ward lower momenta, because the downstream plasma ve- 
locity U2 decreases as the postshock gas cools. The combi- 
nation of the shift of r<,g„ toward lower momenta and the 
decrease of the particles in the Maxwellian tail due to the 
gas cooling leads to the self-regulation of the injection rate 
at a quite stable value. According to the plot of g{p) at 
t = 120 to and t = 240 to, the Maxwell distribution turns 
into a power-law at an almost constant "effective injection 
momentum" which determines the magnitude of the CR 
distribution function g{p) at a stable value (about 1/200 
of the thermal peak). The value of this constant effective 
injection momentum can be translated into the parameter 
/pth ~ 2.3 ( where Pth = Z^mpfceT) defined by 



> 

(D 




[keV] 



Fig. 6. Omni-directional flux vs. proton kinetic energy, for 
t = (dotted), t = lOto (dot-dashed), t = 120 to (dashed), 
and t = 240 to (solid line). These distributions are identical to 
those shown in Fig. 5. For the parameters used see Fig. 4 and 
corresponding text. 



The narrow injection pool also leads to a rather sharp 
transition from the Maxwell distribution to the non- 
thermal part starting shortly above the effective injection 
momentum (see Fig. The canonical result in the test 
particle limit, gip) — p'^ f{p) = p'^p^'^ = constant, for 
a strong shock with s — 3r/(r — 1) = 4 is reproduced 
very well in our simulations. The same energy spectrum 
is shown in Fig. ^in the form of the omni-directional flux 
F{E)dE oc vp'^f{p)dp vs. proton kinetic energy down- 
stream of the shock normalized to t = 0. At energies 
above the injection pool we expect, for the strong shock 
(r ~ 4) simulated here, the result F{E) oc E~°' , with 
a = {{r + 2)/(r — l)}/2 — 1, which is reproduced with 
high accuracy. 

In using the standard cosmic-ray transport equation, 
we have, of course, made use of the diffusion approxima- 
tion, which may introduce an error especially for w ~ M2- 



Cl 



Kang & Jones ( 19951) . But this is somewhat larger than 
what they used (ci = 1.9 — 1.95). 



Using an eigenfunction method. Kirk & Schneider (1989) 
have explicitly calculated the angular distribution of ac- 
celerated particles and accounted for effects of a strong 
anisotropy especially at low particle velocities. They were 
able to calculate the injection efficiency without recourse 
to the diffusion approximation, and found always lower 
efficiencies compared to those in the diffusion approxima- 
tion. Using the initial thermal distribution, we have es- 
timated an effective injection momentum from the peak 
of the distribution function, (7m(p)'''osc(p)- For the shock 
parameters considered here and for e = 0.35 we get an 
effective initial injection velocity of about 6700 km s^^ (in 
the shock frame). For this injection velocity, r = 4 and 
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uq — 5000 km s^^, they estimate a reduction effect of 
« 8%, leaving the diffusion approximation as quite rea- 
sonable even in this regime. 

4-3. Injection and acceleration efficiencies 

To describe the injection efficiency often a parameter ^ 
is used for the fraction of the in-flowing plasma particles 
that are instantaneously accelerated to a fixed injection 
momentum pi„j (e.g. Falle & Giddings 1987; Dorfi 1990] ; 



Jo nes fc Kang 1990 



1993; Berezhko et 



Zank et al. 

al. |1994| ). The injection energy flux / transferred to CRs 
is then given by 



PlUl Pi, 



2m ' 



(18) 



where Ui is the upstream plasma velocity in the shock 
frame, and pi is the upstream density. From the fact that 
the injected energy flux / must be equal to the spatial 
integral of the injection energy loss term S{x), that is, / — 
J S{x) dx and by assuming momentarily a step function 
for the transparency T^^dp) = Hip — Pinj), we get: 



dx 



■ dx , 



ni ui 



(19) 



where ni — pi/m is the upstream number density. This 
is equivalent to the injection parameter used by Kang 



& Jones (1995). The so-deflned injection parameter ^ is, 
however, not an exact measure of the number of particles 
contributing to the population of cosmic rays, because the 
acceleration process cannot be described by shifting par- 
ticles instantaneous from thermal energies to an injection 
momentum Pi„j. Furthermore, ^ depends strongly on the 
chosen injection momentum , which is not a flxed single 
parameter in our numerical simulation. 

A method to measure the injection efficiency with- 
out specifying the injection momentum, is to compare the 
number of particles in the CR part to the number of parti- 
cles swept through the shock. According to our definition 
of the CR population we can write for the CR number 
density 



ncnix,t) = I fc^{p,x,t)d^p 

ifiP,x,t) - f„{p,x,t))d^p . 



(20) 



The fraction of particles that has been swept through the 
shock after the time t, and then injected into the cosmic- 
ray distribution is then given by 



C[t) 



J ricn (x, t)dx 



niUi t 



(21) 



The time development of this injection efficiency is plot- 
ted in Fig. for three valu es of the inverse wave ampli- 
tude e. Recall that Malkov (|199S|) found 0.3 ^ e ^ 0.4. In 
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Fig. 7. Energy efficiency ri(t) and the fraction of cosmic-ray 
particles ^* (t) for three values of the inverse wave-amplitude e 
at a strong, Af — 30, shock. For the parameters used see Fig. 4 
and corresponding text. 



the very beginning of the simulation the injection does de- 
pend strongly on the wave-amplitude, because of the very 
steep dependence of the Maxwell distribution at the injec- 
tion energies. However, as described above, a strong initial 
injection leads to a temperature decrease of the plasma, 
and to a shift of the Maxwell distribution, which balances 
this effect. Therefore at later times the fraction of injected 
particles, does not depend strongly on the initial wave- 
amplitude. At time t/to — 250 (or t = 1.0 ■ lO^s) we get a 
fraction of injected particles of £_* = (1.5 ± 0.4) • 10""^ for 
the interval e = 0.35 ± 0.5. 

To measure the efficiency of the particle acceleration at 
a shock front, we compare the energy ffux in cosmic rays to 
the total energy which is available from the downstream 
plasma ffow. This energy consists of the sum of kinetic 
energy and the gas enthalpy. The fraction of this initial 
energy flux, which is transferred to CRs is given by 



^,u,{t)P,{t) 



^Pdul 



7g 



7g- 



(22) 



where Ud — Us{l — 1/r) is the initial downstream plasma 
velocity in the upstream rest frame. The definition of the 
effic iency r]{t) is similar to the definition from Volk et 
al. (1984). However, Eq. (22) compares the energy flux 
in CRs not only to the kinetic energy flux of the gas, 
but also includes the gas enthalpy flux. We measure the 
CR pressure immediately downstream of the sub-shock. 
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where it will first reach the constant downstream value, 
in case a steady state does exist (see below). The time 
dependent values are averaged over the interval x/xq — 
[— 0.5M2t, • . ■ , 0] in the shock frame to avoid influence of 
small scale modifications of the cosmic-ray pressure and 
plasma velocity on the injection efficiency. When the quan- 
tities U2, Pc and 7c have reached steady-state distribu- 
tions downstream of the sub-shock, rj{t) is also no longer 
time dependent. 

The evolution of the energy efficiency, 77 (i), is plotted 
in Fig. 1^ for three different magnetic-field wave ampli- 
tudes. See Fig. ^ and the description in Sect. ^ for the 
corresponding parameters. The case e = 0.4 corresponds 
to the highest injection efficiency and therefore leads to 
the highest cosmic-ray pressure. To assure a vanishing 
value of the cosmic-ray pressure at the spatial grid bound- 
aries at all times, the calculation for e = 0.4 was done 
on a somewhat larger grid with 60416 uniform zones for 
x/xq = [—300,300]. F or the value e — 0.35, which was 
calculated by Malkov ( 1998 ), we see that about 20% of 
the available energy in this shock is transferred into the 
cosmic-ray population. The acceleration efficiency has, 
however, not reached a real steady state value, but is in- 
creasing with ri(t) oc t" with a ~ 0.1. The acceleration 
efficiency achieved by this time is given by 77 = (18 ± 5)% 
for e = 0.35 ±0.5. Thus a substantial amount of the initial 
energy flux at a shock front can be transferred to a high 
energy part of the distribution, during the relatively short 
time we have simulated here. 

5. Results for weak shocks 

When the initial compression ratio decreases for a weak 
shock, the injection process is influenced in several ways 
by the change in the plasma and magnetic fleld proper- 
ties. To investigate the effects of a lower compression ra- 
tio and lower Mach number on the injection process we 
will consider an example with r = 2.5 and M = 2.24. 
At such a shock, the phase space for which the down- 
stream particles can re-cross the shock to upstream is de- 
creased compared to the strong shock case, because the 
shock velocity in the downstream rest frame U2 = Ug/r 
is inversely proportional to the compression ratio. At the 
same time the plasma is heated less, because the trans- 
formation of kinetic energy to thermal energy depends 
also on the compression ratio; Ak^T oc mpUg(l — l/r^). 
This shifts the downstream Maxwell distribution to lower 
energies, as compared to higher compression, and, there- 
fore, influences strongly the number of particles in the 
momentum range making the potential injection pool. 
On the other hand, at quasi-parallel shocks, the ampli- 
tude of the magnetic fleld wave spectrum Bj^ is amplifled 
downstream by the factor r. For a decreasing compression 
ratio, the downstream plasma becomes more transparent. 
This balances the effects of the phase space and tempera- 
ture changes described above. The initial downstream (in- 
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Fig. 8. Gas density p/po, pressure Pg/Pgo, velocity u/uo, and 
cosmic-ray pressure Pc/Pgo, at times t = (dotted), t = 70to 
(dashed) and t = 140 to (solid line) . The shock Mach number 
is M = 2.24, e = 0.6, and = 100. The initial upstream gas 
pressure is P = 0.12Pgo. 



verse) wave-amplitude e = Bq /B^ was calculated to be in 
the interval e w [0.3, . . . , 0.4] in the limit of strong shocks 



(Malkov |l998t Malkov & Volk |l998| ). An extrapolation to 
weak shocks with r = 2.5 of this interval by multiply- 
ing e with the factor of (4/2.5) gives e « [0.48, . . . ,0.64]. 
However, the calculation of the transparency function was 
based on the assumption of an high amplitude wave spec- 
trum downstream (e <^ 1). With decreasing wave ampli- 
tude the velocity dependence of the transparency function 
changes towards its asymptotic function, deflned by par- 
ticle kinematics without a wave fleld: t{v) = for ti < M2 
and t{v) = 1 for ii > U2- On the other hand, this limit 
may be reached in reality only if the resulting beam from 
downstream to upstream is too weak to produce a mag- 
netic fleld instability. 

As an initial exploration of this behavior, we will 
present here results for the spatial and momentum dis- 
tributions and the energy and particle injection efflciency 
for an inverse magnetic flelds amplitude parameter e in the 
range e = [0.4, . . . , 0.7] . We have included the value e = 0.4 
to compare the results directly to the strong shock case. 
This can demonstrate the principal effect of weaker shocks 
on the injection process. The resulting injection efficien- 
cies and shock modifications for all values of e shown here 
should be considered as lower limits for the weak shock 
with r = 2.5 {M = 2.24) as described above. 

The physical scales are specified as follows: to ~ 1.11 • 
10^ s, xo = 3.33 • 10^3 cm, uq = 3000kms"\ po/wp = 
0.03 cm-3, = 4.52 • lO^^ergcm-^. We use C = 100 
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Fig. 9. Phase-space density g = f vs. proton momentum 
immediately downstream of the sub-shock. Also shown is the 
transparency function Tcsc- Both functions are presented for 
t = Q (dotted), t = lQto (dot-dashed), t = 70to (dashed), and 
t = 140 to (solid line). For the parameters used see Fig. 8 and 
corresponding text. 



for the simulations presented here, and a magnetic field 
of -B = 3/iG. The initial values for the M = 2.24 case 
are pup = Po, Wup = -uo, and Pg,up = 0.12Pgo in the 
upstream region, while pd = 2.5po, U2 — — 0.4mo, and 
-Pg.d = 0.72Pgo in the downstream. We have used 44032 
uniform grid zones for x/xo ~ [—170, 130], with the shock 
initially at x = 0, and 128 uniform grid zones in log(p) 
for log(p) = [—3.0, 0]. The corresponding Mach number is 
M = 2.24. 

Figure^ shows the normalized gas density p{x), gas 
pressure P^{x), plasma velocity u{x) and the cosmic- 
ray pressure Pc{x) over the spatial length x, for different 
times. Because the resulting non-thermal spectrum pro- 
duced as a result of the injection and particle acceleration 
is steeper than in the strong shock case, the pressure Pc in 
this distribution remains small compared to the gas pres- 
sure at all times. As a result, the shock is modified only 
slightly. Also the temperature of the downstream plasma 
remains almost constant. Furthermore, because the en- 
ergy density in non-thermal particles is not an increasing 
function in time, the shock modification can reach a steady 
state earlier, as compared to the strong shock case. In fact, 
at time t = 140to, shown in Fig. ^, the pressure Pg, Pc, 
the velocity u and the density p immediately downstream 
has reached almost a steady state. 

The downstream momentum distribution in Fig. |9| 
shows clearly the steeper spectrum of the non-thermal 
part, which asymptotes to the standard result g{jp) oc 
p-^+'i with s = 3r/{r — 1) = 5 for r = 2.5. It can be 
seen also, that the thermal part of the distribution is 
not as much modified as in the strong shock case (com- 
pare Fig. H) . Because the modification of the transparency 
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Fig. 10. Energy efficiency ri{t) and the fraction of cosmic-ray 
particles ^*{t) for four values of the inverse wave-amplitude 
e at a weak shock. For the parameters used see Fig. 8 and 
corresponding text. 



function over time depends only on changes in the down- 
stream plasma velocity, it remains essentially unchanged. 

The energy efficiency 'ri{t), as defined in Eq. (22), is 
lower roughly by a factor of two compared to the strong 
shock case, because of the steeper non-thermal spectrum 
and the resulting energy density (compare Fig. |^ and 
Fig. |l^. Our results for the wave amplitude, 0.5 < e < 
0.7, give the injection efficiency, ^* = (2.5 ± 0.7) • 10^^ 
at time t = lAQto = 1.55 • lO^s, where the time evolution 
can be considered as almost a steady state. The number 
of particles, which are in the non-thermal part is compa- 
rable to the strong shock considered above at this time. In 
addition, we point out that the application of the above 
described injection model to weak shocks is an extrap- 
olation, and we believe would yield lower limits on the 
injection efficiency. 

6. Conclusions 

We have developed a numerical method to include self- 
consistently the injection of the supra-thermal particles 
into the cosmic-ray population at quasi-parallel shocks 
according to the analytic solution of Malkov ( 1998 ). 
Toward this end, we have adopted the "transparency 
function" t^sc{v,U2) which expresses the probability that 
supra-thermal particles at a given velocity can leak up- 
stream through the magnetic waves, based on non-linear 
particle interactions with self-generated waves. We have 
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incorporated the transparency function into the exist- 
ing numerical code which solves the cosmic-ray transport 
equation along with the gas dynamics equations. In or- 
der to investigate the interaction of high energy parti- 
cles, accelerated by the Fermi process, with the under- 
lying plasma flow without using a free parameter for the 
injection efficiency, we have applied our code with the 
new injection scheme to both strong {M = 30) and weak 
(M = 2.24) parallel shocks. 

The main conclusions from the simulation results are 
as follows: 

1. The injection process is regulated by the overlap of 
the population of supra-thermal particles in the in- 
jection pool and the function dT^^^{p,t)/dp. As be- 
ing in the high energy tail of the Maxwell velocity 
distribution, the population in the injection pool de- 
pends strongly on the gas temperature and the particle 
momentum. The function dT^,^{p,t) /dp behaves like a 
delta-function defined near a narrow injection pool. As 
the postshock gas cools due to high initial injection, the 
Maxwell distribution shifts to lower momenta. But the 
transparency function also shifts to lower momenta, as 
well, due to its dependence on the postshock flow ve- 
locity. As a result, the injection rate reaches and stays 
at a stable value after a quick initial adjustment, and 
also depends only weakly on the initial conditions. This 
self-regulated injection may imply a broad application 
of our simulation methods. 

2. The fraction of the background particles that are accel- 
erated to form the non-thermal part of the distribution 
turns out to be in the range 1.2 • 10~^ ^ 1-9 • 10~^ 
for the range of initial wave-amplitudes 0.3 < e < 0.4 
at a M = 30 shock. For a M = 2.24 shock, a slightly 
higher injection is achieved at ^* — (2.5 ± 0.7) • 10"^, 
but this could be a lower limit. Such values for the par- 
ticle injection efficiency have been used as a parame- 
ter for spherically expanding SNRs by several authors 
(Dorfi |199G| ; Jones fc K ang |l992t Berezhko et al. |1995| ; 



Berezhko & Volk 2000). These values are well above 
the "critical injection rate" of ?7crit ^ lO""* above which 
spherical shocks of this Mach number are CR domi- 
nated according to Berezhko et al. ( |1995| ). 

3. Due to computational limitations of using a Bohm 
type diffusion model, we have integrated our mod- 
els until the maximum momentum reaches about 
(Pmax/"ipc) ^ 1. For the M = 30 shock model, the 
energy flux in the total CR distribution was about 
18% ±5% of the energy flux in the thermal plasma and 
shocks didn't become CR dominated and smoothed 
completely by the end of our simulations. For the 
M = 2.24 shock model, the acceleration efficiency is 
lower by a factor of two compared to the high Mach 
shock because of the smaller velocity jump across the 
shock. 

4. Just above the injection pool, the distribution func- 
tion changes sharply from a Maxwell distribution to 



an approximate power-law whose index is close to the 
test-particle slope. We estimated this critical momen- 
tum as Pi„j ~ (2.2 - 2.3) • where = 
This determines the number of particles in the injec- 
tion pool by /(Pinj) oc exp{—pf^^./2mpk^T). For strong 
shocks this translates into a distribution function at 
injection energies of g{p,^^) (1/100- l/200)g(pth)- 
While the weak shock model of M — 2.24 reaches a 
steady-state, the strong shock model of M — 30 has not 
reached a steady-state by the end of our simulation. We 
expect for the strong shock that the CR pressure continues 
to increase and the shock becomes CR dominated, lead- 
ing to the greater total velocity jump and more efficient 
acceleration. In realistic shocks such as SNRs, however, 
escaping particles due to non-planar geometry or lack of 
scattering at high momentum are likely to become impor- 
tant. To resolve this non-linear evolution, much longer 
physical time scales have to be simulated, until CRs reach 
energies where escape is likely to be important. The key 
problem here is the range in configuration and momen- 
tum space that has to be computed. Our method uses a 
grid with uniform cells in configuration space, chosen fine 
enough to capture the evolution of g{x,p) at near-thermal 
momenta where the diffusion coefficient is proportional to 
p^ (Bohm diffusion). This leads to a computationally ex- 
tremely expensive calculation, especially because the grid 
has to be large enough to contain the diffusion length scale 
of the highest momentum CRs. The problem can be solved 
on a much larger time scale by using an adaptive mesh re- 
finement (AM R) cod e with the shock tracking techniques 
(Kang & Jones 1999 ) . In the near future we plan to incor- 
porate the injection model presented here into the pow- 
erful shock tracking AMR-code, to calculate the evolu- 
tion of the phase-space distribution of the plasma during 
different phases of SNRs. This would allow us to inves- 
tigate with a plasma-physical based injection model how 
the slowly growing cosmic-ray pressure at a strong shock 
eventually modifies the shock structure. A strong modifi- 
cation will cause the velocity jump across the subshock to 
decrease and the distribution function of the suprather- 
mal particles to steepens. This might have further back 
reaction on the injection efficiency. Also the CR distri- 
bution will deviate from a simple power-law. For a cal- 
culation up to the highest energy CRs, also the spher- 
ical geometry of a SNR should be taken into account. 
Such an approach could lead to a consistent calculation 
of the complete phase-space distribution at quasi-parallel 
shocks, and should be a promising step towards a calcu- 
lation of the overall efficiency of SNRs in producing CRs 
during their evolution. 

For oblique shocks, the injection efficiencies calculated 
here for a parallel shock should define an upper limit, be- 
cause the statistical probability of a particle to cross the 
shock from downstream to upstream decreases with the in- 
tersection velocity of magnetic field and shock front . This 
kinematical effect was investigated by Baring et al. (1993) 
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with the use of Monte-Carlo simulations. However, in the 
model we have incorporated here, the injection is already 
suppressed strongly (compared to the purely kinematical 
model) by the reduced transparency of the plasma due to 
the high amplitude Alfven waves. We point out, that for 
oblique shocks, the filtering due to Alfven waves may be 
reduced due to the decreased downstream amplification of 
the wave amplitude. This would allow lower energy par- 
ticles to be injected, and the kinematical effect could be 
partly balanced. As a result, we speculate that the depen- 
dence on the obliquity might be significantly weaker than 
calculated by Baring et al. ( 1993 ). Resolution of that im- 
portant question must await more complete understanding 
of the injection physics. 

In summary, we have shown that the process of parti- 
cle acceleration under consideration of a plasma physical 
injection model underlies a rather effective self-regulation. 
Apart from the direct particle-wave interaction described 
by the injection model itself, also the energetic feedback of 
the energy transfer between thermal plasma and cosmic- 
rays keeps the fraction of particles in the non-thermal dis- 
tribution at roughly 10^"^ of the particles swept through 
the shock. These self-regulation mechanisms lead to a 
quite stable injection efficiency, which depends weakly on 
the initial conditions. 
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